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In a recent contribution [Phys. Rev. B 81, 165104 (2010)] fermionic Projected Entangled-Pair 
States (PEPS) were used to approximate the ground state of free and interacting spinless fermion 
models, as well as the t-J model. This paper revisits these three models in the presence of an 
additional next-nearest hopping amplitude in the Hamiltonian. First we explain how to account for 
next-nearest neighbor Hamiltonian terms in the context of fermionic PEPS algorithms based on sim- 
ulating time evolution. Then we present benchmark calculations for the three models of fermions, 
and compare our results against analytical, mean-field, and variational Monte Carlo results, re- 
spectively. Consistent with previous computations restricted to nearest-neighbor Hamiltonians, we 
systematically obtain more accurate (or better converged) results for gapped phases than for gapless 
ones. 
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I. INTRODUCTION 

Several recent papers have proposed and explored the 
use of tensor networks to simulate fermionic lattice mod- 
els in two spatial dimensions)^ including algorithms 
based on the Multi-scale Entanglement Renormaliza- 
tion Ansatf^l (MERA) and Projected Entangled-Pair 
States EHUD The fundamental ingredient, common to all 
approaches, is to incorporate fermionic statistics directly 
into the tensor network by regarding the tensors as lin- 
ear maps of anticommuting degrees of freedom, an idea 
recently generalized to anyonic statistics! 17 * 18 * The main 
goal of fermionic tensor network methods is to address 
strongly correlated fermionic models, which suffer from 
the negative sign problem in Quantum Monte CarloJ^ 

Fermionic PEPS were first proposed in Ref. [2] and 
have been discussed in several other papers In Ref. [71 
we provided a detailed account of how to adapt ex- 
isting bosonic PEPS algorithms to the fermionic case, 
and we us ed th e fermionic version of the infinite PEPS 
algorithm d 13 * 16 ! to obtain benchmark results for three 
models with nearest-neighbor Hamiltonian in an infinite 
square lattice: (i) a model of free spinless fermions with 
a pairing potential, (ii) a model of interacting spinless 
fermions with a nearest-neighbor repulsion, and (Hi) the 
well-known t-J model. In the first case, a comparison of 
the numerical results with the exact solution showed that 
fermionic PEPS could reproduce ground-state energies 
and short-range correlators satisfactorily, with a larger 
degree of accuracy in gapped phases than in gapless ones. 
For the model of interacting spinless fermions, PEPS 
yielded significantly lower variational energies than ob- 
tained by mean-field theory (restricted Hartree-Fock the- 
ory) , which enabled to determine the phase diagram more 
accurately. For the t-J model, PEPS energies are com- 
parable (or even better in some cases) than variational 
Monte Carlo based on Gutzwiller-projected ansatz wave 
functions. 

As with any new approach, systematic benchmarking 



of fermionic PEPS algorithms is important in order to es- 
tablish their range of applicability. The results of Ref. [71 
while limited to three specific models, were a first step 
in this direction. A key question to be addressed is how 
good a fermionic PEPS is in practice, as a variational 
ansatz, at approximately representing the ground state 
of fermionic models. Of course, the precise answer to this 
question will depend on the specific model under consid- 
eration. However, insight on how fermionic PEPS meth- 
ods generally perform in certain circumstances, e.g. in a 
given gapped phase, may be obtained from models where 
an exact solution or previous numerical results by other 
methods are already available. Such insight is essential 
in order to subsequently assess the validity of fermionic 
PEPS results obtained in more relevant (and challenging) 
scenarios, such as in exploring the ground state phase di- 
agram of the t-J model, which was addressed in Ref. [6j 
or of the Hubbard model. 

Thus, one of the main goals of this paper is to fur- 
ther benchmark the performance of fermionic PEPS algo- 
rithms, by considering more complex models than those 
addressed in Ref. 7. Specifically, here we will consider 
the effect of adding nearest-neighbor hopping terms to 
the three models of Ref. |7[ Recall that in many cases of 
interest it is desirable to consider a model where fermion 
particles can hop between nearest-neighbor sites (with 
amplitude t) as well as next-nearest neighbor sites (with 
amplitude t'). For example, in effective models of high- 
T c superconductors (cuprates), it is estimated that the 
ratio |t7t| is of the order of 0.1 - 0.3P*mJ a fi nite t > 
can have several important effects on th e system. For 
instance band-structure calculation d 21 * 22 * and experimen- 
tal analysis 2 ^ suggest that the highest T c strongly de- 
pends on t' ft. Previous studies of the hole-doped t- 
t'-J model revealed that a finite t' < can suppress 
magnetic orde r 2 4 " 26 ! and enhance or suppress pairing 
correlation d 26 * 27 * (depending on the doping). It was a lso 
shown that t r influences the formation of stripes pi!29] 

In order to explore how well fermionic PEPS can repro- 
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duce such ground states, we need to extend the fermionic 
PEPS algorithm of Ref. 7J based on simulating imaginary 
time evolution, to the case where the evolution is gener- 
ated by a Hamiltonian that also contains next-nearest 
neighbor terms. Thus a second main goal of this paper 
is to explain how this is accomplished. The simulation of 
frustrated spin models with next-nearest neighbor terms 
with PEPS by imaginary time evolution was considered 
in Ref. [30l In contrast to Ref. [30j here we will explain 
how to generalize the so-called simple update scheme for 
time evolution to the case of next-nearest neighbor terms 
(in addition to accounting for the fermionic character of 
the PEPS). 

The paper is organized as follows: In Sec. [II] we first 
explain how to update the PEPS during an imaginary 
time evolution in the presence of next-nearest neighbor 
Hamiltonian terms. Then we discuss how to evaluate the 
expectation value of a next-nearest neighbor operator, as 
required e.g. in order to compute the energy of a PEPS 
in the case of nearest neighbor hopping t' . In Sec. Ill we 



present a series of benchmark results for extensions, con- 
taining next-nearest neighbor hopping terms, of the three 
models addressed previously in Ref. [71 namely (i) an ex- 
actly solvable model of free spinless fermions with t' ^ 0, 
(ii) the t-t'-V modelPand (in) the t-t'-J modelPThe 
accuracy and/or apparent convergence (as a function of 
bond dimension D) of ground state properties in these 
models are comparable to tho se p reviously obtained for 
the case t' = 0. Finally, Sec. IV summarizes our find- 



ings and conclusions. For completeness, Appendix [A] pro- 
vides details on the two different corner-transfer-matrix 
(CTM) schemes used in order to evaluate expectation 
values from an infinite PEPS. 



II. METHOD 

As in Ref. 13 for spin systems, we use an infinite PEPS 
with bond dimension D to approximate the ground state 
of a Hamiltonian defined on an infinite square lattice, by 
simulating an evolution in imaginary time starting from 
some (random) initial state. 

The evolution itself is first approximated, through a 
Trotter-Suzuki decomposition, by a sequence of two-site 
gatesP^ After applying each two-site gate, the affected 
bond in the PEPS has to be truncated, so that the 
evolved state is again represented by a PEPS with the 
same bond dimension D. This truncation implies choos- 
ing a D-dimensional subspace in the vector space asso- 
ciated to the bond index. In Ref. we distinguished 
between two different truncation or update schemes. 
The first consists of choosing the subspace that best 
supports the wave function. This requires taking the 
whole PEPS wave-function into account during the up- 
date, i.e. the environment has to be computed at ev- 
ery step in the imaginary time evolution.^ Alternatively, 
following Refs. [15] one can update the PEPS as in the 
time-evolving block decimation (TEBD) method in one 



dimension^ by means of a singular value decomposition 
(SVD). In this second, simpler option, referred to as sim- 
ple update, instead of considering the full environment, 
only local weights attached to the bonds of the PEPS 
are taken into account. In one dimensional systems with 
open boundary conditions, this choice of update is op- 
timal, since the full environment can be encoded in the 
local weights. In two dimensional systems, however, the 
simple update is no longer optimal (a better chose of 
truncated space can be obtained by considering the whole 
environment) but it has a significantly lower computation 
cost as a function of bond dimension D, which allows to 
consider larger bond dimensions (in a suboptimal way) 
and potentially obtain more accurate results with the 
same computational cost. 

Here we will use the simple update. In this section 
we first discuss the additional steps needed to apply the 
simple update of Ref. [7] in the case of a two-site gate 
acting on next-nearest neighbor sites. Then we will also 
explain how to evaluate the expectation value of a next- 
nearest neighbor operator. 




FIG. 1: (Color online) (a) Infinite PEPS with a four-site unit 
cell with tensors A, B, C, and D. (b) Representation of the 
same PEPS with diagonal matrices on the bonds, which 
is used for the simple update, (c)-(f) Relation between the 
tensors in (a) and (b). 



A. Simple update for next-nearest neighbor terms 

In order to perform the simple update the fermionic 
PEPS in Fig. [TJa) is recast into the form shown in 
Fig. [ljb) , where the diagonal matrices live on the 
bonds and the tensors T q live on the sites of the network. 
As already explained in Appendix B in Ref. the sim- 
plified update for nearest-neighbor links consists of the 
three steps summarized in Fig. [2] 

The next-nearest neighbor update is performed in a 
very similar way, with the difference that three PEPS 
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(a) 




FIG. 2: (Color online) Simple update of a nearest- neighbor 
link, (a) A gate g is applied to a link between tensors Fa 
and Ts. Contracting the diagram, including all adjacent ma- 
trices Afc, yields tensor G. (b) Singular value decomposition 
of tensor G leads to tensors Fa, Fb and diagonal matrix Ai. 
(c) The updated tensors F' A , F' B are obtained by multiplying 
the corresponding unaltered inverse diagonal matrices A^ 1 to 
Fa and as shown. The new diagonal matrix X[ corre- 
sponds to the D largest diagonal entries (singular values) of 
Ai. 



tensors and two diagonal matrices are updated at the 
same time, through two consecutive singular value de- 
compositions. Figure [3] illustrates the update for the link 
between tensors Fa and Fd, via the tensor Fb (including 
all adjacent diagonal matrices A&). Similarly one could 
also perform the update for the same link involving the 
tensor Fc instead of Fb as shown in Fig. [4^a). In prac- 
tice we apply the square root of the gate to both com- 
binations of tensors, in order to make the update more 
symmetric. In principle it is conceivable that also the or- 
der in which the singular value decompositions are made 
plays a role. However, in the cases studied we have not 
found a significant difference when changing the order. 

Figures Qc)-(d) show the relevant diagrams for the 
update of the other diagonal link in the tensor network. 
As usual, crossings in the network have to be replaced 
by swap tensors (black diamonds) in order to account for 
the fermionic anticommutation rules. In total there are 8 
different diagonal links for the 2x2 unit cell, where each 
link is updated in two different ways (via two different 
intermediate tensors). 

We conclude this section with three remarks. Firstly, 
we note that the complexity of the update can be reduced 
by splitting off the parts of the tensors involved in the up- 
date by a singular value decomposition (see e.g. Fig. 32 in 
Ref.[7j). Secondly, the same update may of course be used 
also for bosonic and spin systems, with the simplification 
that crossings do not need to be taken into account, since 
bosonic and spin operators commute. Finally, note that 
an update for next-to-next nearest neighbor interaction 



(a) 




FIG. 3: (Color online) Simple update of a next-nearest neigh- 
bor link, (a) A gate g is applied to a link between tensors 
Fa and Fd, via the tensor Fb- Contracting the diagram, 
including all adjacent matrices Xk, yields tensor G. (b) Sin- 
gular value decomposition of tensor G leads to tensors Fa, G', 
and Ai. (c) Tensor O" includes the weights Ai, obtained by 
keeping the D largest diagonal entries of Ai. Singular value 
decomposition of tensor G ; leads to tensors Fb, Fd, and A4. 
(d) The updated tensors F' A , F' B , and F' D are obtained from 
Fa, Fb, and as shown. The new diagonal matrix A 4 cor- 
responds to the D largest diagonal entries of A4. 



can be implemented in a very similar way, since it also 
involves three PEPS tensors (arranged on a line). This 
case, not further considered in the paper, requires a larger 
unit cell of size 3x3. 



B. Computing expectation values of next-nearest 
neighbor terms 

In order to compute the expectation value of a next- 
nearest neighbor operator, one proceeds in a similar way 
as explained in Sec. Ill B of Ref. 7 for nearest neigh- 
bor operators. First, one has to compute the environ- 
ment S^-cdI for the tensors A, £?,C, D, which accounts 
for the infinite lattice surrounding the 2x2 unit cell 
formed by these tensors. In Ref. [71 this was done with 
the directional corner transfer matrix (CTM) method.^ 
Besides this scheme, in the present work we also use 
another variant of the CTM approach, as discussed in 
Appendix [Aj The CTM algorithm yields the four cor- 
ner tensors Ci, 62,63,(74 and the eight edge tensors 
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SVD 



FIG. 4: (Color online) The remaining relevant diagrams of 
the next-nearest neighbor updates. The updated tensors are 
obtained in a similar way as explained in Fig. [3] 



Tzi, T r i, T w2 , T d2 , T Z3 , T r3 , T n4 , T dA shown in Fig. [5| which 
altogether constitute the environment £^cd\ 

Next, one connects the four tensors A, £?, C, D together 
with their complex conjugates to the environment, and 
joins the physical legs to the operator accordingly, as 
exemplified in Fig. [5] for a next-nearest neighbor opera- 
tor acting between A and D. Since the wave function 
encoded by the iPEPS is not normalized, the value ob- 
tained by contracting the tensor network in Fig. [5] has 
to be divided by the norm of the iPEPS, which is sim- 
ply obtained by replacing the two-site operator o by the 
identity operator in Fig. [5] 

Evaluating a two-site operator o linking tensors £?, C 
can be done in a similar way, simply by reconnecting the 
legs of the operator o (highlighted in green in Fig. [5| 
to the physical legs of B and C accordingly. Finally, 
the expectation value of the remaining six next-nearest 
terms can be obtained analogously by first generating the 
environments for the plaquettes [ p C ] , [ ^ ] , [ ^ ] . 



are chosen randomly. In some cases we observe that the 
state resulting from the time evolution depends on the 
initial condition (seed), which might be due to the local 
character of the simple update (cf. Sec. II A). We there- 
fore typically run several (of the order of 10) simulations 
with different seeds, and pick the state with lowest en- 
ergy. In some cases a good choice is to initialize a PEPS 
from a previously converged PEPS with smaller D, but 
this does not always lead to the state with lowest en- 
ergy. A converged PEPS for a certain Hamiltonian usu- 
ally provides a good initial condition for simulations with 
slightly different Hamiltonian parameters, provided both 
states are in the same phase. 

Evaluating the expectation value of observables re- 
quires computing an environment, which in turn requires 
introducing a secon d bond dimension x associated to ad- 
ditional truncations! 7 * 13 * 1 ^ In all present simulations the 
bond dimension x °f the environment has been chosen 
to be sufficiently large so that the expectation values of 
local observables do not significantly change when fur- 
ther increasing x- Typical values are x = 36 for D = 2, 
X = 48 for D = 4, and x = 64 for D = 6 and D = 8. 



A. Free fermions including a pairing potential 

We first consider a model of free spinless fermions given 
by the Hamiltonian 



III. BENCHMARK RESULTS. 

In this section we provide benchmark results for three 
fermionic models on an infinite square lattice with near- 
est neighbor and next-nearest neighbor hopping terms, 
with amplitudes t and t' . Each model is an extension of 
an analogous model with only nearest neighbor hopping 
term, that is with t' = 0, previously addressed in Ref. [Tj 

As initial condition, the tensors of the infinite PEPS 



+ * E^+#.c.]-2A^4q, (1) 

(«» i 

with (ij) and ((ij)) denoting the sum over nearest and 
next-nearest neighbor pairs, respectively. In the follow- 
ing we fix the hopping amplitude t to 1 and the pairing 
potential 7 to 1, and consider a chemical potential A in 
the range [1,4]. For t' = the model reduces to the one 
studied in Sec. IV of Ref. [7j and is gapped for A > 2 and 
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critical for A < 2.^ For t' ^ 0, the location of the tran- 
sition Ac between gapped and gapless phases depends 
on t'. 
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FIG. 6: (Color online) Relative error of the ground state en- 
ergy of the free spinless fermion model as a function of 
A, for different values of D and t' . Full and open symbols 
correspond to critical and gapped phases, respectively. 



Figure [6] shows the relative error of the ground state 
energy as a function of A for different D and t' . As in the 
case for t' = the energies are improved upon increasing 
D, and the accuracy is higher in the gapped phase (open 
symbols) than in the critical phase (full symbols). In 
general the error is larger than in the t' = case (cf. 
Ref.[7|), but still of the order of 1(T 3 (1(T 5 ) in the critical 
(gapped) phase for D = 6. 

Figure [7] shows the two-point correlation function 



C(r) 



(2) 



as a function of distance r between the two sites (in x- 
direction). The numerical results are seen to approach 
the exact values with increasing D, with correlations at 
short distances being better reproduced than correlations 
at long distances (see middle panels). Also in this case, 
the accuracy is better in the gapped phase (7 = 1,A = 
3.5, t' = 0.4) than in the critical phase (7 = 1, A = 2, t' — 
0.6). 



B. The t-t'-V model 

Next we study the t-t'-V spinless fermion model and 
compare the infinite PEPS results with the mean-field 
studies from Ref. |3TJ based on Hartree-Fock (HF) theory 
restricted to states invariant under translations by two 
sites. The Hamiltonian reads 

(ij) (ij) 

- E^^+^ c -]-^E^, (3) 
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FIG. 7: (Color online) Upper panels: Correlation function 
C(r) — (cjci+r) as a function of distance r (in x-direction) 
in the gapless (left) and gapped (right) phase of the free- 
fermion model (IT]). Middle panels: Absolute value of C(r) in 
semi- logarithmic scale. Lower panels: The difference between 
the simulation result C(r,D) and the exact result C ex (r) for 
different values of D. 



with V being the nearest-neighbor interaction strength, 
and fx the chemical potential. As an example we study 
the transition between a metal phase at low electron den- 
sity n to a charge- density- wave (CDW) phase at half fill- 
ing (n = 0.5), for fixed parameters t = 1, t 1 = —0.4, and 
V = 2. [A similar study was done in Ref. 7 in the case 
of t' = 0.] 

The HF study predicts a first order phase transition be- 
tween the metallic phase for densities < n < 0.274(1) 
and a CDW phase, which is thermodynamically stable 
upon doping in the range 0.368(1) < n < 0.5, as il- 
lustrated in Fig. [8^a). The region in 0.274(1) < n < 
0.368(1) is unstable, corresponding to phase separation 
(PS) between the two states. The study was done for 
systems of size 100 x 100. 

Figure [9] shows a comparison of the energies in the 
two phases, obtained with HF and infinite PEPS. As ex- 
plained in Ref. [7] for the case t' = 0, the crossing of the 
two energies is obtained by starting from a state deep 
in the metal (CDW) phase and then increase (decrease) 
H across the transition. Similarly as for t' = 0, PEPS 
energies in the gapped CDW phase (for n = 0.5) do not 
differ significantly from the HF energies. However, with 
increasing D = 4, 6, 8, PEPS energies do differ signifi- 
cantly from HF energies in the metal phase close to the 
transition. This produces a shift of transition point /1* 
to larger values of /i, corresponding to a larger value of 
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the density n, as shown in Fig. [8[b). Also, in contrast 
with the HF prediction, which predicts a stable doped 
CDW phase, we do not find a stable doped CDW phase 
for D > 4, i.e states in the doped CDW phase obtained 
by our PEPS computation exhibit a higher energy than 
states in the metal phase or in the CDW phase at half fill- 
ing. Thus, we only find a stable CDW phase at exactly 
half filling (for the model parameters under considera- 
tion). 
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FIG. 8: (Color online) (a) Mean- field phase diagram for fixed 
parameters V = 2 and t' — —0.4 as a function of particle den- 
sity n, obtained by Hartree-Fock (HF) restricted to states in- 
variant under translations by two sit esP^ The charge-density- 
wave (CDW) phase is separated from the metal by a region 
of phase separation (PS). The CDW phase is stable upon 
doping, b) Phase diagram obtained with iPEPS for different 
bond dimensions D for the same parameters as in a) . A stable 
CDW phase is only found at half filling. The phase boundary 
to the PS region is shifted towards higher values of n with 
increasing D. 



The amount of entanglement in the gapped CDW 
phase is relatively low, as can be seen in the lower panel 
in Fig. [9j showing the entanglement entropy of a 2 x 2 
block in the system. Therefore, a PEPS with small bond 
dimension is already sufficient in order to obtain an ac- 
curate description of the ground state. In the metal 
phase, however, the entanglement entropy is considerably 
higher. The energies have not yet converged as a func- 
tion of the bond dimension D for D = 8, and thus further 
corrections to the energy can be expected for larger val- 
ues of D. This phase appears to be particularly difficult 
to represent by the PEPS. This is not surprising. In 
the presence of a ID Fermi surface the entanglement en- 
tropy exhibits a logarithmic multiplicative correcti on to 
the area law, as shown in the case of free fermions! 33 * 34 * 
A PEPS representation, however, can only reproduce an 
strict area law of the entanglement entropy,^ and there- 
fore it is unclear that the ansatz can offer an accurate 
approximation of the ground state of a metallic phase. 
Nevertheless, our results also show that, with increasing 
bond dimension D, the PEPS can still be used to obtain 
a systematic improvement over HF results. Since the 
phase boundary does not change much when comparing 
the D = 6 with the D = 8 simulations, it seems likely 
that the D = 8 result is already close to the exact one. 

We conclude this section with two remarks. Firstly, we 
observed states in the metal phase close to the transition 
that exhibit a slight density modulation between sublat- 
tices A and B in the lattice, similarly as in the CDW 
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FIG. 9: (Color online) Upper panel: Energy per site of the 
t-t'-V model g as a function of chemical potential fi for 
V — 2 and t' = —0.4, obtained with (restricted) Hartree-Fock 
(HF) 31 and infinite PEPS. The first order phase transition be- 
tween the metal phase and the charge- density wave (CDW) 
phase occurs at a value fi* where the two corresponding en- 
ergies cross. Middle panel: Particle density n as a function 
of chemical potential in the two phases. At the first order 
phase transition point /i*, n jumps from a certain value n* 
in the metal phase to n — 0.5 in the CDW phase. For den- 
sities in between n* and n = 0.5 the system exhibits phase 
separation. In contrast to the HF study, infinite PEPS sim- 
ulations do not yield a stable doped CDW phase for D > 4. 
Lower panel: Entanglement entropy of a 2 x 2 block in the 
two phases, illustrating that there is substantially more en- 
tanglement in the metal phase than in the CDW phase. 



phase. However, this modulation becomes weaker with 
increasing D, which strongly suggests that this symmetry 
breaking of translation invariance is a numerical artifact 
due to small D rather than a real physical feature. 

Secondly, it is conceivable that the t-t'-V model ex- 
hibits also CDW phases other than the checkerboard- 
ordered phase at half filling, i.e. with a period larger 
than 2. Indeed, in some simulations we observed states 
where the density of e.g. sublattice A oscillates as a func- 
tion of the CTM steps, which could be an indication for 
a CDW phase with a period larger than 2. However, in 
order to represent such states accurately by a PEPS, ei- 
ther a rather large bond dimension or a larger unit cell 
than the ones employed in this work would be required. 
We therefore point out that the final phase diagram is 
likely to be different than the one presented in Fig [8] 
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t-t'-J model 



IV. DISCUSSION AND CONCLUSION 



Finally, we present benchmark results for the energy 
of the t-t'-J model, given by the Hamiltonian 

H = - * E [ g U> + H.c] + J^2(SiS j - ^%)(4) 

(ij)o- (ij) 

with a = {T, 4-} the spin index, hi = J2 a cl a Ci a the elec- 
tron density and Si the spin 1/2 operator on site z, and 

Qcr — Q<j(l ^ia^*^)" 

Here we compare our results of the energy with 
variational Monte Carlo (VMC) and fixed-node Monte 
Carlo (FNMC) results from Ref. [26j which are based 
on Gut z wilier-projected ansatz wave functions including 
spin and density Jastrow factors. 

Figure [To] shows the energy as a function of particle 
density n for J/t = 0.4 and t' jt — —0.2. Similarly as 
in the case t' = of Ref. 7, the results for D < 8 have 
a higher energy than VMC, but for D = 8 the energies 
are comparable or even lower than VMC. Thus, the ad- 
ditional t' does not seem to change the accuracy of the 
ansatz significantly (compared to Fig. 26 in Ref. 7). The 
FNMC results, however, are still considerably lower than 
the D = 8 results. Note that the Monte Carlo ener- 
gies are for a finite lattice size (with 98 and 162 sites), 
and that the energy increases with increasing system size. 
Thus, the FNMC energy in the thermodynamic limit is 
slightly higher than shown in the plot. Finally we point 
out that for D = 8 the maximal dimension x we used is 
64. It is conceivable that the energies still change slightly 
when increasing x further, but we do not expect a signif- 
icant change. 
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FIG. 10: (Color online) Energy per site (with the chemical 
potential term subtracted) as a function of particle density 
n of the t-t'-J model, with J/t = 0.4 and t' /t = -0.2. The 
best iPEPS results for D = 8 lie in between the VMC and 
the FNMC results. 



In this paper we have explained how to extend the 
fermionic PEPS algorithm for infinite lattices presented 
in Ref. 7 to models including next-nearest neighbor terms 
in the Hamiltonian. This opens the possibility to study 
a much larger variety of models, which are relevant e.g 
for high-temperature superconductivity.^ 

The benchmark results in Sec. |III[ involving three dif- 
ferent models, indicate that the accuracy of the ground 
state energy is comparable to the case where only nearest- 
neighbor terms are present in the Hamiltonian. The 
present results are compatible with, and provide addi- 
tional evidence in favor of, the conclusions of Ref. [71 
where it was indicated that a fermionic PEPS seems to 
offer a more accurate description of fermionic for gapped 
phases than for critical phases, among which two types 
(I and II) need to be distinguished. 

Gapped phases -Out simulations produced high accu- 
racies (or better convergence as a function of bond di- 
mension D) for ground state energies for gapped systems, 
such as the free spinless fermion model ffl) for large A and 
the interacting spinless fermion model (pi at half filling. 
The corresponding ground states exhibit a comparatively 
small amount of entanglement and can therefore be ap- 
proximated accurately by a PEPS with small bond di- 
mension D. 

Gapless phases of type /.-Gapless systems with a finite 
number of zero modes in their spectrum (that is, with- 
out a ID Fermi surface) are in general more entangled 
than gapped systems, but they are believed to still obey 
the area law of the entanglement entropy, which a PEPS 
is known to be able to reproduce.^ Therefore, a PEPS 
is still expected to offer an accurate description of the 
ground state, provided certain bond dimension D, in gen- 
eral larger than in the gapped case, is used. This category 
of states includes, for the free spinless fermion model ([!]), 
the gapless p-wave paired phase corresponding to small 
A; and, for the t-t'-J model, both the antiferromagnetic 
phase at half-filling and the (expected) d-wave paired 
phase in the doped case. A remarkable achievement 
of fermionic PEPS simulations is that they yield better 
or comparable energies than the usual Gutzwiller pro- 
jected ansatz wave functions for the doped t-t'-J model 
for D = 8. However, a larger bond dimension D would 
be needed to attempt to match the energies of state-of- 
the-art fixed-node Monte Carlo (FNMC)P 1 

We note that even if the energy in such gapless phases 
is obtained with a few digits of accuracy (of the order 
of 0.1% for the model with D = 6), it is still un- 
clear whether the PEPS reproduces other relevant prop- 
erties of the ground state accurately. Here we found that 
the correlation function ([2| for the free fermion model 
([I]) is satisfactorily reproduced for short-range distances. 
However, there are other known cases, e.g. the Heisen- 
berg antiferromagnet model, where the accuracy of the 
order parameter is two orders of magnitude worse than 
the accuracy of the energy for a D = 5 PEPS on an 
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infinite latticed Therefore, the reliability of PEPS re- 
sults must be carefully checked on a case by case ba- 
sis. Nevertheless, and taking into account that there 
are no exact methods available to address systems of 
strongly correlated fermions, we believe that fermionic 
PEPS and, more generally, fermionic tensor networks, of- 
fer a useful approach to such systems that complements 
other approaches such as fixed-node Monte Carlo J™ dia- 
grammatic Monte Carlof^ cluster dynamical mean-field 
theory^ or Gaussian Monte Carlo!™ 

Gapless phases of type 77-Gapless systems with a ID 
Fermi surface, such as the metal phase of the interacting 
spinless fermion model ([3|, are known to display loga- 
rithmic multiplicative corrections to the area law of the 
entanglement entropy. This logarithmic violation can 
not be reproduced with a PEPS, and it is therefore un- 
clear that fermionic PEPS methods will be able to ac- 
curately describe the ground state of such massively en- 
tangled phases. Nonetheless, our results for the metal 
phase of the interacting spinless fermion model ([3| show 
that, even in this case, PEPS with increasing values of 
the bond dimension D can be used to obtain system- 
atic improvements on mean-field energies, which in turn 
question the validity of the mean-field phase diagram. 

It might well be, however, that a proper characteri- 
zation of the ground state of gapless phases with a ID 
Fermi surface is simply beyond the reach of fermionic 
PEPS. In this case other techniques, such as one of the 
many methods which work particularly well in Fermi- 
liquid type phases at weak coupling (e.g. Refs. 37 39j) , or 
a specialized tensor network approach™ should be used 
instead. 

We conclude by noticing that a larger bond dimen- 
sion D, and therefore more accurate PEPS results, may 
be within reach in subsequent studies. This larger val- 
ues of D could be accessed e.g. by using more computer 
resources (possibly in a parallel architecture) , by exploit- 
ing internal the symmetries of the the fermionic models^ 
(e.g. particle conservation) and/or by employing Monte 
Carlo sampling techniques. 42 We also point out that with 
the same values of D = 2 — 8 used in this paper, more ac- 
curate results may be obtained by employing the standard 
update instead of the simple update in the simulations, 
which, however, comes with a larger computational cost 
(cf. Ref.EJ). 
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Appendix A: Variant of the corner transfer matrix 
method 

The evaluation of the expectation value of a local ob- 
servable from a PEPS requires the computation of the 
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FIG. 11: (Color online) Left-right renormalization step of the 
anisotropic CTMRG method, (a) Tensor network of the re- 
duced tensors a,6,c, and d embedded in the environment. Cut- 
ting the lines as marked in the figure leads to tensors Qi and 
Q2 shown in (b). (b) A singular value decomposition of the 
tensors Qi and Q2 is performed, and only the x largest singu- 
lar values are kept, (c) Tensor network obtained by inserting 
two extra rows of tensors in the middle of the tensor network 
of (a), and by introducing an approximate resolution of the 
identity I ~ U\U\ in (each of the the two copies of) cut 1, 
and similarly for cut 2. [One could also use W\ instead of Ui]. 
The cuts 3 and 4 define tensors Q3 and Q4, and a singular 
value decomposition of these tensors is performed, similarly 
as in (b), yielding tensors U3 and U4. (d) Renormalized cor- 
ner and edge tensors. Note that the top and bottom edge 
tensors have to be swapped after the update, i.e. T[ x — T r i, 

T»/ rri rri I rri rri/ rri 

rl — J- 11, -LIZ — l-r3, lr3 = J- 13- 



so-called environment ] 7 ^ 13 ^ which accounts locally for 
the ground state wave function on the rest of the system, 
as explained in Sec. |IIB| In an infinite system, corner 
transfer matrix (CTM) methods, originally introduced 
by Baxter J™ can be used to approximately compute this 
environment. In this case, the environment is approxi- 
mated by four corner tensors Ci, C2, C3 and C4 and eight 
edge tensors (or half-row transfer matrices) T/i, T ri , T n2 , 
Td2, T t3 , T r3 , T n4 , and T d4 , as shown in Fig. [5] 

In the present work we applied two different CTM 
schemes: The first is the directional CTM methodp^ al- 
ready used in our previous workP the second scheme 
essentially corresponds to the CTM renormalization 
group (CTMRG) method from Ref. HI adapted to the 
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anisotropic case and to a 2 x 2 unit cell, with subse- 
quent coarse-graining moves in horizontal and vertical 
direction as in the directional CTM method. The main 
difference is that in the former scheme the four corners in 
the lattice are renormalized individually (see Refs. I7|16l 
for details) , whereas in the CTMRG approach the full en- 
vironment is taken into account in each renormalization 



step. Figure 11 illustrates the left-right coarse-graining 
move, i.e. where the system size is increased by two lat- 
tice sites in the horizontal direction. In a similar way a 
top-bottom move is performed, which increases the sys- 
tem by two lattice sites in the vertical direction. These 
two moves are iterated until convergence is reached. In- 
cluding more tensors in each renormalization step helps 
to better determine the relevant subspace to be kept dur- 
ing truncation, since more information about the system 
is taken into account. 1 The disadvantage is that this 
scheme is more strongly limited by the machine precision 
1 A similar idea is pursued in the second renormalization group 
method.^ 



of the computer. This can be understood by considering 
the isotropic case: In the directional CTM of Ref. [16] 
the spectrum of a single corner matrix C is computed, 
whereas in the present scheme, which involves multiply- 
ing all four corner matrices, yields the fourth power of the 
same spectrum, with singular values (in this case, eigen- 
values) expanding many more orders of magnitudes, and 
therefore more vulnerable to errors due to finite machine 
precision. 

We observed that the CTMRG scheme converges bet- 
ter in the case of highly-entangled systems, as for ex- 
ample in the metal phase of the model (|3| close to the 
phase transition. In gapped systems both methods seem 
to converge equally well, with the directional CTM yield- 
ing slightly better accuracies than the present scheme. 
We used the former to address the gapped phase of the 
model ([!]), and the latter in all other simulations. 
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